import numpy as np
import differential as p1
import matplotlib.pyplot as plt

def resolution_eq(y0,t0,N,h,meth):
    """ Resolution des equations de Malthus et Verhuslt """
    gamma = 2.9
    k = 20e9

    Malthus = lambda y,t : y
    Verhulst = lambda y,t : gamma*y*(1-(y/k))
    X = p1.meth_N_step(y0,t0,N,h,Malthus,meth)
    Y = p1.meth_N_step(y0,t0,N,h,Verhulst,meth)
    plt1, = plt.plot(X)
    plt2, = plt.plot(Y)
    plt.legend([plt1,plt2],["Malthus","Verhulst"])
    plt.show()

#resolution_eq( 20e9 , 0. , 100 , 1/float(20) , p1.step_rk4 )
#resolution_eq( 40e9 , 0. , 100 , 1/float(20) , p1.step_rk4 )

def resolution_Lotka_Volterra(y0,t0,N,h,meth):
    """ Resolution du systeme et affichage de solutions autour du point de depart ainsi que des variations """
    a = b = c = d = 1
    
    F = lambda y,t: np.array([ y[0]*(a-b*y[1]) , y[1]*(c*y[0]-d) ])

    for i in range(0,5):
        res = p1.meth_N_step(y0+i*0.1, t0, N, h, F, meth)
        X = []
        Y = []
        for j in res:
            X.append(j[0])
            Y.append(j[1])

        plt3 = plt.plot(X , Y)
    plt.title("Solution autour du point de depart")
    plt.show()
    
    plt1, = plt.plot(X)
    plt2, = plt.plot(Y)
    plt.legend([plt1,plt2],["Proies","Predateurs"])
    plt.title("Variations")
    plt.show()
    return X,Y

#resolution_Lotka_Volterra(np.array([0.5,0.05]), 0, 1000, 1/float(100), p1.step_rk4)
#resolution_Lotka_Volterra(np.array([0.5,0.5]), 0, 1000, 1/float(100), p1.step_rk4)
#resolution_Lotka_Volterra(np.array([0.5,5]), 0, 1000, 1/float(100), p1.step_rk4)

def calcul_periode(X):
    """ Recherche les deux premiers pics de X et retourne la difference entre leurs indices """
    n = len(X)
    i = 1
    while(  not( (X[i-1]<= X[i]) and (X[i] >= X[i+1]) ) and (i < n-1)  ):
        i += 1
    p1 = i

    i+=1

    while(  not( (X[i-1]<= X[i]) and (X[i] >= X[i+1]) ) and (i < n-1)  ):
        i += 1
    p2 = i
    return p2-p1


resolution_eq( 6e9 , 0. , 100 , 1/float(20) , p1.step_rk4 )
print calcul_periode(resolution_Lotka_Volterra(np.array([0.5,0.5]),1,1000,1/float(100),p1.step_rk4)[0])
